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ABSTRACT 

Observations of quasars at z > 6 suggest the presence of black holes with a few times 10® M©. Nu¬ 
merous models have been proposed to explain their existence including the direct collapse which provides 
massive seeds of 10^ Mq. The isothermal direct collapse requires a strong Lyman-Werner flux to quench 
H 2 formation in massive primordial halos. In this study, we explore the impact of trace amounts of metals 
and dust enrichment. We perform three dimensional cosmological simulations for two halos of > 10® Mq 
withZ/Zo = 10 10“^ illuminated by an intense Lyman Werner flux of J 21 = 10^. Our results show that 

initially the collapse proceeds isothermally with T ~ 8000 K but dust cooling becomes effective at densi¬ 
ties of 10* - 10*^ cm“* and brings the gas temperature down to a few 100-1000 K for Z/Zq > 10“®. No 
gravitationally bound clumps are found in Z/Zq < 10“^ cases by the end of our simulations in contrast to 
the case with Z/Zq = 10“"^. Large inflow rates of > 0.1 Mo/yr are observed for Z/Zq < 10“^ similar to a 
zero-metallicity case while for Z/Zq = 10 the inflow rate starts to decline earlier due to the dust cooling 
and fragmentation. For given large inflow rates a central star of ~ 10"^ Mq may form for Z/Zq < 10^^. 

Subject headings: methods: numerical - cosmology; theory - early Universe - galaxies: formation 


1. Introduction 


The discovery of high redshift quasars at z > 6 
reveals the presence of supermassive bl ack holes 

(SMBHs) of about a few billion solar masses ^fin et alJ _ 

2006: ^illott et al j2010l: Mortlock et alj201 u Venemans et al 

fb 


2OI3I: IWu et all I2OI5I: IVenemans etZI l2015h. The 

formation of SMBHs a few hundred million years 
after the Big Bang presents a challenge for our un¬ 
derstanding of early structure formation. The path¬ 
ways for their assembly include the stellar mass 
black holes forming from the collapse of Pop III 


stars I 


(I Abel et alJl2002t Yoshida et'dll2008: Latif et al. 


2013a Hirano et al.l l2014l) . black hole seeds from 
the collapse of a dense nucl ear cluster either due 


to th e relativistic instability (iBaumgarte & Shaniro 
1999), via stellar dynamical proc e ss dOmukai et al. 


2008; Devecchi & Volonteril 20091 : Katz et alT 20151: 


Yaiima & Khochfail l2015h ~ even the core collapse 
of a cluster consisting of stellar mass black holes 


(iDavies et alJl201 It iLupi et al.ll20l4) or a direct col¬ 
lapse of a protogalactic ga s cloud dLoe b & Rasio 
1994: Bromtn & LoebI 2003 : iBegelman et al. 20061: 


Spaans_&_Sil^l20oJ Latifet^ 2ol3a|ldyFOTarae^^ 


20141: Inavoshi et al.l2014l: Choi et al.l2015t Shlosman et al. 
20161). Details of these models are discussed in dedi 


cated reviews dVolonteri 2010l : Volonteri & Bellovarv 


2012l:lHaimanH2012h. 

The masses of seed black holes from the above 
mentioned scenarios range from a few hundred to 
thousand solar masses except the direct collapse model 
which provides massive seeds of about 10^ - 10^ Mq. 
The lighter seeds require continuous accretion close 
to the Eddington limit in order to reach a billion so¬ 
lar masses in a few hundred million years. There¬ 
fore, massive seeds forming via the direct collapse are 
favoured for the assembly of high redshift quasars. 
The formation of so-called direct collapse black 
holes (DCBHs) requires large inflow rates of about 
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> 0.1 Mrri/vr dBeselman 2010l: Hosokawaet al. 2013: 

2007: Wise et alJ 20121 Safranek-Shraderet al. 2014: 

Schleicher et al.l 201. 3l Sakurai et al.l 20L5|). These 

Bovino et al. 

2014 Ritter etalj 2015). The critical 


inflow rates can be obtained either via dynamical 
processes such a s the ’bars-i n-bars’ instability (see 
Begelman et al. ( 2006h and Begelman & ShlosmanI 


(120091) 1 or thermodynamically due to a large speed 
of sound (M ~ c^/G ~ 0.1 Mo/yr(T/8000 
where Cs is the thermal sound speed). Large inflow 
rates can be more easily achieved in an isothermal 
direct collapse which necessitates the presence of a 
strong Lyman Werner (LW) flux to quench the for- 
mation of molecular hydrogen in met a l free halos 
(Omucai 200 it Schleicher etal. 2010t Shang et al 


2010t 


Latif et al 


2014) . The critical value of the LW 


flux (Jcrit) above which isothermal collapse occurs 
depends on the spectra of the stars (Omu^ 200l|_ 
Shang etal. 2010t Latif et al. 20141: Sugimura et al. 


2014: Agarwal & K iochfaill 20L5 : Regan et al. 2014bl: 

Latif et alT 201.5 ). Latif etal. (l20L5h computed the 

strength of Jcrit for realistic Pop II spectra from three 
dimensional cosmological simulations and found that 
it corresponds to a value of J 21 of a few times lO'* 
where J 21 = 1 implies a flux of 10“^' erg/cm^/s/Hz/sr 
below the Lyman limit. Another constraint for the 
isothermal direct collapse is that the halos should be 
metal free as trace amounts of metals and dust can lead 


to fra gmentation and star formation (lOmukai et al. 

2nn8h. 


Numerical simulations show that large inflow rates 
of about 0.1-1 Mo/yr are available in the massive pri¬ 
mordial halos of 10^ - 10* Mq illuminated by a strong 
LW flux at z=15 and massive objects of 10^ Mn can be 
formed within ab ou t 100,000 years (|W ise_eL^ 200^ 
Latifet^ 20]Aa||^ReganeL^ 20]4ayBeceme^^ 


2015 : Inavoshi et al.l l20l4 Van Borm et alJ 20141) . 


However, the value of Jcrit (> 10"^ in terms of J 21 ) 
required for the formation of isothermal DCBHs 
is much higher than the background UV flux and 


ties of star forming galaxies (^Diikstra et alJ 2008t 

Affarwal et alJ 20121 Diikstra et al. l20 

4 Visbal et al. 

2014: Latif et al. 2015; Habouzit et al. 

2 OI 5 I). Suner- 

nova winds from nearby galaxies ca 
halos with metals and make the sites 
mal direct collapse even more rare 

n also pollute 
for an isother- 

(.Johnson et al. 


2013h . In the context of second generation star for¬ 


mation (Pop II), it has been found that above a critical 
value of metallically (Zcut = 3 x IO^^^Zq) fragmen¬ 
tation becomes inevitable and leads to the formation 
of Pop II stars ( Omukai et alJl200-5 ; Glover & .lappsenl 


value of the metallicity is further reduced by about two 
orders of magnitude in the pres e nce of trace amounts 


of dust ( Schneider et al. 2003 : Omukai et al. 20051 


20081) and th e H? formation on dust grains enhances 
the cooling dOrnukai 2001 : Cazaux & Spaans 20091: 
Latif et al. 2012h . Numerical simulations show that 
low mass star formation seems plausible due to the 
dust cooling occurring at high densities in minihalos 
with metallicities as low as a few times 10~^ solar 
j Tsuribe_&Omu^ 2006[ 2008; Dopcke et al. 2013 
Smith et alJl20L5h 


Omukai et al. ( 20081) explored the thermodynami¬ 


cal properties of gas irradiated by a strong LW flux 
above Jcrit and enriched by trace amounts of metals 
and dust. They found from a one-zone model that 
in such conditions dust cooling becomes effective for 
Z/Zo > 5 X10 ® and proposed that a dense stellar clus¬ 
ter can form under these conditions which may later 
collapse into a black hole of 100-1000 Mq. 

In this study, we explore the impact of dust and 
metal line cooling in massive primordial halos illumi¬ 
nated by a strong background LW flux and polluted by 
trace amounts of metals and dust. To accomplish this 
goal, we perform high resolution three dimensional 
cosmological simulations for two halos of a few times 
10^ Mq by turning on a LW flux of strength 10^ in units 
of J 21 at z=30 and presume that they are pre-enriched 
with Z/Zq = 10“® - 10“"^. Our results suggest that 
dust cooling occurs even for Z/Zq = 10 ® at densities 
above 10'*^ cm“* but still the conditions are suitable 
for the formation of massive objects. This work has 
important implications for assessing the feasibility of 
black hole formation models. 

This article is organised as follows. In section 2, we 
describe numerical methods and chemical model em¬ 
ployed in this work. We present our results in section 
3 and discuss its implications for black hole formation 
in section 4. Summary of our main results and conclu¬ 
sions are discussed in section 5. 

2. Numerical Methods 


We have performed simulati ons using the ope n 
source code Enzo version 2.4 ( Brvan etaP 2014 ). 
Enzo is an adaptive mesh refinement, parallel, grid 
based cosmological simulation code which has been 
extensively used to perform high resolution simula¬ 
tions. It uses the message passing interface (MPI) 
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to achieve scalability and portability on various plat¬ 
forms. The piece-wise parabolic method (PPM) is 
employed to solve hydrodynamics, the particle-mesh 
technique (PM) is used to solve dark matter dynamics 
and the multi-grid Poisson solver for gravity. 


2.1. Simulation setup 


Our simulations are started with cosmological ini¬ 
tial conditions at z=100 generated with the ’’inits” 
package available with Enzo. The computational do¬ 
main has a size of 1 Mpc/h in comoving units and pe¬ 
riodic boundary conditions are adopted both for grav¬ 
ity and hydrodynamics. We use parameters from the 
WMAP 7 years d ata to generate initial conditions 
( .larosik et alJl2011 ). Our computational box is cen¬ 
tered on the most massive halo forming at z=15 se¬ 
lected from uniform grid DM only simulations of res¬ 
olution 128^ particles. We rerun the simulations with 
a top grid resolution of 128^ grid cells (same number 
of DM particles) and add two nested refinement lev¬ 
els each with a resolution of 128^ grid cells and 128^ 
DM particles. In all, we employ 5767168 particles to 
solve the DM dynamics which provide us an effective 
DM resolution of about 600 Mq. We smooth DM par¬ 
ticles at about a 2 pc scale to avoid spur ious numerical 
artefacts, see iLatif & Volonteril (12015h for a detailed 
discussion. We further add 28 dynamical refinement 
levels during the course of the simulations in the cen¬ 
tral 62 kpc comoving region of the box. This approach 
enables us to resolve the gravitational collapse down 
to scales of about 0.1 AU. We employ a fixed Jeans 
resolution of 32 cells throughout the simulations to re¬ 
solve turbulent eddies and to fulfil the Truelove crite¬ 
rion which requi res a resolution of at-least four cells 
per Jeans length (Truelove et al.ll 19971 : iFederrath et al. 


2011 : Latif et al. 2013bl) . Our dynamical refinement 
criterion further includes the baryonic over density and 
particle mass resolution. Grid cells are marked for re¬ 
finement if gas density exceeds four times the cosmic 
mean. Similarly, cells are flagged for refinement if DM 
density is above 0.0625 times where poM is the 

dark matter density, r = 2 is the refinement factor, ( is 
the refinement level, and a — -0.3 makes the refine¬ 
ment super-Lagrangian. 

We select two halos of 5.6 x 10^ Mq and 3.25 x 
10^ Mq. The properties of these halos such as their 
spins, collapse redshifts and Jcrit are listed in table 1 of 
Latif & Volonteri ( 2015 ) as the halo A and the halo C, 
respectively. We presume that the halos are enriched 
with trace amount of metals and dust by supernova 


winds from a nearby star forming galaxy. The lat¬ 
ter also provides a strong LW flux above the critical 
threshold required for the quenching of H 2 formation. 
The prime aim of this study is to explore the impact 
of dust cooling in the presence of an intense LW flux. 
We therefore select J 21 = 10^ for which halo collapses 
isothermally in the absence of dust cooling. At den¬ 
sities above 10“*® g/cm-* where dust cooling becomes 
important the gas cloud is already self-shielded against 
external radiation and therefore the value of Jait does 
not depend on the metal content of the halo. This 
scenario is particularly relevant for the halos which 
get polluted while being irradiated by a strong LW 
flux, and we specifically explore here the cases with 
Z/Zo = 10“®, 10“® and 10“^. 


Due to the computational constraints, we stop our 
simulations when a peak density o f 3 x 10“'* g/cm^ 
is reached. Omukai et al.l ( 2008h have shown that 
dust cooling becomes optically thick at densities 
above ~ 10“*** g/cm^ for the metallicities explored 
here, and thereafter collapse is expected to pro¬ 
ceed adiabatically. The transition to the adiabatic 
regime during the gravitational collapse of primor¬ 
dial gas has been explored in hydrodynamical sim¬ 
ulations ( InavpshietaL 2014; Van Borm et al. 2014 


Becerra et alJ 2015 : Latif et al. 2016h . These studies 
found that the adiabatic phase stabilises the grav¬ 
itational collapse on small scales, breaks its self¬ 
similarity and a pr otostar begins t o form at this stage. 
For dust cooling, lOmukai et al.l (120081) have shown 
that the transition to the adiabatic phase occurs at ear¬ 
lier stages of the collapse (at lower densities compared 
to the primordial gas). They also considered the sput¬ 
tering of dust grains in their model finding that the dust 
temperature remains below the vaporisation tempera¬ 
ture in optically thin regime and once adiabatic phase 
begins, it exceeds the evaporation temperature but gas 
cloud is already optically thick. Therefore, thermal 
evolution remains adiabatic similar to the metal free 
gas. Furth er details of our model are described in 


Omukai et al.l (120081) . 


In this study, we explore three cases with metallic¬ 
ities of Z/Zq = 10“®, 10“® and 10“^* for each halo and 
turn on a strong LW fl ux of strength J 21 = 10®, well 
above the Jcrit found in lLatif et al.l (l2015l) . Further de¬ 
tails of our chemical network as well as heating and 
cooling processes are described in the section below. 


3 































































2.2. Chemical model 


We employ the KROME package (IGrassi et al. 


20141) to solve the chemical and thermal evolution of 
the gas in cosmological simulations. The rate equa¬ 
tions of H, H+, H , He, He+, He++, H 2 , H+, e are 
solved to study their evolution during the gravitational 
collapse. We employ a uniform background LW flux 
of strength J 21 = 10^ emitted from a Pop II galaxy 
which can be mi micked by a rad i ation t emperature 


of 2 X 10^ K (see ISugimura et al.l (120141) : iLatif et al 


(l20I5h ). Our results do not depend on the choice of 
the radiation spectrum as long as the flux is well above 
the critical thre shold to dissociate molecular hydro¬ 


gen (d erived by lOmukai et al.l (l2008h : ISugimura et al 


(2014) for the one-zone model and by iLatif et al 


(12015h in 3D simulations). The list of chemical reac¬ 


tions i s given in table 1 in the appendix of ILatif et al. 


(I 2 OI 5 I) . Our model includes H 2 formation, H 2 and 
photo-dissociation, photo-detachment, collisional 
induced emission and H 2 collisional dissociation. In 
addition to this, we add a reacti on for H 2 forma tion on 
dust grain surfaces taken from Omukair ( 2000l) . Vari¬ 
ous cooling and heating mechanisms such as cooling 
due to the collisional excitation, collisional ionisa¬ 
tion, radiative recombination, Bremsstrahlung radi¬ 
ation and chemical cooling/heating from three-body 
reactions are included in the chemical model. We also 
employ the H 2 self- s hieldi ng fitting formula given in 
Wolcott-Green et al. ( 2011 ). 

In addition to the cooling and heating processes 
mentioned above, we include cooling due to the 
fine structure lines of CII and OI, cooling and heat¬ 
ing by dust grains which is relevant for Z/Z© = 
10“^, 10“^ and 10 Our treatment for metal and dust 
cooling/heating proc esses comes from the model of 
Omukai et al. ( 2008h . here we summarise their main 
features. We do not explicitly include the chemical re¬ 
actions involving metals but instead assume that most 
of the oxygen and carbon are in OI and CII form due 
to the presence of a strong LW flux. Carbon has a 
lower ionisation energy (11.26 eV) than hydrogen and 
therefore is in CII form while OI remains in the atomic 
state due to the charge exchange with hydrogen atoms 
in a neutral medium. The dust grain composition and 
size distribution is assumed to be similar to the so¬ 
lar neighbourhood and the amount is scaled with the 
metallicity of the gas cloud. We assume a dust to gas 
mass ratio of 0.01 Z/Z©. Dust grains in the supernova 
ejecta are more efficient in cooling and H 2 forma¬ 
tion due to their smaller size and larger surface area 


( Schneider et al. 2003 : Omukai et alJ 20081) . The dust 
is assumed to be in thermal equilibrium and its tem¬ 
perature (Tdust) is evaluated separately from the gas 
temperature at each density and temperature. The fol¬ 
lowing e nergy balance equat ion is solved to compute 
the Tdust ( Omukai et al.ll200^ 


4;r 


^ difst) dv — A.gas^dust 

J" Ka,vBviTcMB)dv. 


( 1 ) 


+47: 


Here Agas^dust is the net gas cooling rate per unit mass 
due to dust grain collisions, Ka^y is the dust absorption 
opacity an d By is the black body spectr um. Agas^dust 
is given as ( Hollenbach & McKeell 19791) 


1/2 


( T \ ^ 

1 ^) ^ 

( —15K\ 

(T-Tdus,)(ZIZ^) 


where n is the gas number density and Tdust is the dust 
temperature. We have set temperature floor at Tcmb- 
Lurther details about the dust proper ties can be found 
in section 2.2 of Omukai et al. (2008). 


3. Results 

We here present our main results obtained both 
from a one-zone model and three dimensional cosmo¬ 
logical simulations. 

3.1. One-zone model 

To test the chemical model presented in the pre¬ 
vious section, we performed a one-zone test with a 
LW flux of strength J 21 = 10^ (well above the Jcrit 
from one-zone model with T^d = 2 x lO'* K) and 
varied the metallicity from Z/Z© = 10“® - 10^“*. We 
took an initial temperature of 160 K, a gas density of 
~ 10“^^ g/cm^ and species abundances of 2 X 10“® 
for e“ and 2 x 10 for H 2 . In Ligure[T] we show the 
temperature and the abundances of H 2 , H“ and e“ as 
a function of density. In the presence of a strong LW 
flux, H 2 gets photo-dissociated, its abundance remains 
below 10“® and is unable to cool the gas. The temper¬ 
ature increases up to about 10“* K where atomic line 
cooling becomes effective and cools the gas down to 
8000 K up to the densities of 10“^^ g/cm^. At higher 
densities for Z/Z© > 10“^, dust cooling comes into 
play and sharply cools the gas down to a few hundred 
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K. For the case of Z/Zq = 10 H 2 formation on 
dust grains becomes effective and its abundances get 
boosted very quickly. The small increase in the tem¬ 
perature at a density of 10“'^ g/cm^ is due to the chem¬ 
ical heating produced by the formation of H 2 on grain 
surfaces. Dust cooling continues at high densities up 
to 10“'° g/cm^ and after that the cloud becomes opti¬ 
cally thick to dust self-absorption. Beyond this point 
the thermal evolution proceeds adiabatically and the 
inclusion of dust cooling breaks the a self-similarly 
of an isothermal collapse. A protostar is expected to 
form at this stage with a central adiabatic core and is 
expected to grow rapidly by accretion from the proto- 
stellar envelope. 

For Z/Zq = 10“^, the dust cooling becomes effec¬ 
tive at densities of 10“*^ g/cm^, about two orders of 
magnitude higher than for Z/Zq = 10“"'case. The tem¬ 
perature sharply declines down to a few hundred K and 
the H 2 fraction gets boosted. The cloud becomes opti¬ 
cally thick to dust cooling similar to the Z/Zq = 10“"' 
case at 10“'° g/cm^ and the thermal evolution follows 
the adiabatic equation of state. At an even lower metal- 
licity of Z/Zq = 10“°, the dust cooling becomes neg¬ 
ligible and the thermal evolution follows the atomic 
cooling track. The H 2 fraction is increased up to about 
a few times 10“"'. The abundances of e“ and H“ re¬ 
main very low and are similar at high densities for 
Z/Zq = 10“° and Z/Zq = 10“"' cases. The degree 
of ionisation is higher for Z/Zq = 10“°. The impact 
of metal line cooling such as CII and OI is negligible 
for such metalli cities. Overall , our r esults are in good 
agreement with lOmukai et al.l OOOSh . 


3.2. 3D simulations 

In total, we have performed six cosmological sim¬ 
ulations for two halos of 5.6 x 10^ Mq (HI) and 
3.25 X 10^ Mq (H2) for metallicities of Z/Zq = 
10“"', 10“° and 10 °. The collapse redshifts of the 
halos are 10.4 and 10.8, respectively, and independent 
of the dust content. We also performed one additional 
simulation with zero metallicity for the sake of com¬ 
parison for HI. The gravitational collapse of the halos 
was followed from 1 Mpc down to sub-AU scales with 
the adaptive mesh refinement approach. We assumed 
that the halos are irradiated by a strong LW flux of 
str ength 10° in term s of J 21 , well above the JcUt found 


3.2.1. Thermodynamical and dynamical properties 

The density-temperature phase diagrams of both 
halos for Z/Zq = 10“"', 10 ° and 10 ° are shown in 
figure |2] In the presence of a strong LW flux which 
dissociates H 2 molecules, the halos are unable to col¬ 
lapse via molecular hydrogen cooling and continue to 
grow via merging and accretion until they reach the 
atomic cooling limit. Both halos have recently gone 
through a major mer ger and the merger histor y of these 
halos is described i mLatif & Volonteril ( 2015h . At den¬ 
sities of about > 10“^"' g/cm°, atomic line cooling be¬ 
comes effective and cools the gas to about 8000 K. The 
collapse proceeds isothermally mainly via atomic hy¬ 
drogen cooling up to densities of 10“'° g/cm° and the 
role of metal line cooling remains negligible. At den¬ 
sities > 10“'° g/cm°, dust cooling comes into play 
and brings the gas temperature down to a few hun¬ 
dred K depending on the amount of dust. For Z/Zq = 
10 °, the collapse remains isothermal up to densities 
of about 10“'^ g/cm°, at higher densities the dust cools 
the gas down to a temperature of ~1000 K for halo 2 
while both warm and cold phases coexist for halo 1 at 
densities > 10“'^ g/cm° due to the local variation in 
the gas density and collapse velocity. In the interme¬ 
diate case of Z/Zq = I0“°, the dust cooling becomes 
important earlier at densities of ~ 10“'"' g/cm° for both 
halos and decreases the gas temperature down to about 
700 K. 

For Z/Zq = 10“"', the dust grain cooling becomes 
strong at densities of ~ 10“'° g/cm° and sharply brings 
the gas temperature down to a few hundred K. The 
temperature stalls at about 1000 K at densities between 
~ 10“'"' - 10“'^ g/cm° due to the chemical heating 
from the H 2 formation on grain surfaces and again the 
dust cooling takes over at higher densities. These re¬ 
sults are in agreement with the one-zone test presented 


in the previous section and also with lOmukai et al 


Latif et al.l (120151) . and the LW flux as well as the 


fixed metallicities of the above mentioned amount are 
turned on at z=30. 


(I2008h . Small differences in the thermal evolution are 
observed for both halos, the cooling at high densities 
is generally more effective for halo 2 in comparison 
with halo 1. We show the spherically averaged pro¬ 
files of density, H 2 fraction and temperature in figure|3] 
The collapse remains isothermal down to the scales of 
10,000 AU irrespective of the metallicity explored in 
this work due to the strong LW flux. Below this scale 
dust cooling becomes important and remains confined 
to the central 10 AU for Z/Zq = 10 °, the inner ~100 
AU for Z/Zq = 10 ° and extends up to 3000 AU for 
Z/Zq = 10“"'. The central temperature is about 300, 
700 and 1000 K for Z/Zq = 10“"' - 10“°, respectively. 
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Similarly, the H 2 fraction remains low at scales above 
10^AU and exceeds 10^^ (above which H 2 cooling be¬ 
comes important) within the central part of the halo 
due to its formation on dust grains. For Z/Zq = 10 
the H 2 abundance gets boosted within the central few 
1000 AU while for Z/Zq = 10“^ only in the central few 
hundred AU and for Z/Zq = 10“^ it remains confined 
to the central 10 AU. 

The maximum density reached in our simulations is 
about 3 X 10“" g/cm^ and the density increases with 
R“' ^ above ~ 100 AU, close to the expected profile 
from an isothermal collapse (i.e. The den¬ 

sity profile is shallower for Z/Zq = 10 in the cen¬ 
tral 1000 AU due to efficient dust cooling while for 
the other cases it almost follows the isothermal pro¬ 
file. Small bumps in the profile indicate the formation 
of additional clumps. The turbulent velocity is about 
20 km/s close to the viral radius and increases towards 
smaller scales up to about 30 km/s. This increase in 
the turbulent velocity is due to the infall of gas to¬ 
wards the center and is higher for dust cooling cases 
compared to the metal free case. To quantify whether 
a disk-like structure may form, we computed the ra¬ 
tio of the velocity dispersion (i.e., where 

Cs is the sound speed and Vtmb is the turbulent velocity) 
to the rotational velocity, corresponding to the vertical 
turbulent support within the disk. The latter provides 
an estimate for the ratio of vertical scale height to the 
disk radius. We found that this ratio remains about 1 
for all cases except for halo 2 with Z/Zq = 10“^ where 
Vro,/cr ~ 0.3. This suggests that the rotational veloc¬ 
ity is almost comparable to the velocity dispersion and 
a very thick disk-like structure with H/R of about 1 is 
formed in all cases. To further estimate the rotational 
support against gravity, we calculated the ratio of ro¬ 
tational velocity to the Keplerian velocity (Vrot/VKep)- 
The ratio of Vrot/VKep is shown in figure|4]and is about 
0.6 for all cases. It gets enhanced between 60-1000 
AU for dust cooling cases. A peak in Vrot/Vxep for 
Z/Zq = 10“^ case comes from the density structure in¬ 
side the halo. Overall, the estimates of Vrot/VKep sug¬ 
gest that the rotational support gets enhanced earlier 
for higher metallicity cases. However, enhanced rota¬ 
tional support can delay the runaway collapse but does 
not stop it. Moreover, the rotational support is only im¬ 
portant in the envelope of the halo and decreases down 
within the central core. 


3.2.2. Mass inflow rates 


Our estimates for the mass inflow rates (AnR^pVrad) 
are shown in figure |3] Mass inflow rates of 0.1 - 
1 MQ/yr are observed for Z/Z 0 < 10“^ down to about 
100 AU while for Z/Zq = 10“"^ the mass inflow rate 
starts to decline around 10^ AU for halo 2. For com¬ 
parison, we also show a zero-metallicity case where 
cooling is mainly due to atomic line cooling. We found 
that apart from small differences, the inflow rates for 
Z/Zo < 10“^ almost follow the zero-metallicity case. 
Below 100 AU the mass inflow rates decline in all 
cases. It comes from the fact that during the runaway 
collapse a core-envelope structure develops, where the 
core has flat density and the envelope has a ~ R“^ 
profile, with a core-envelope boundary around 10-100 
AU. The radial velocity is high inside the envelope, de¬ 
creases down in the core and consequently a decline in 
the instantaneous inflow rate occurs towards the cen¬ 
tre. In the central core thermal pressure balances the 
gravity as the core is not yet gravitationally unstable 
which results in a decline of the inflow rate. How¬ 
ever, this may change at later stages when the central 
core further collapses. Fragmentation observed in the 
Z/Zq = 10“"* cases may also partly contribute to the 
decline in the inflow rates. We also note that mass in¬ 
flow rate scales with sound speed cube (M oc cl) and 
therefore decreases with the gas temperature. The dust 
cooling reduces the gas temperature which results in 
lower inflow rates. The accreting matter inside the en¬ 
velope will collapse into the protostar at later times and 
high inflow rates are expected to be maintained. The 
increase in the mass inflow rates was observed in our 


previo us simulations for a metal free case (iLatif et al. 


2013ah as more mass collapsed into the core. We ex¬ 


pect a similar trend for other cases until the feedback 
from the central star blows away the gas or rotational 
support inhibits the mass accretion. 


A similar trend is observed in the radial infall ve¬ 
locity which peaks around 15 km/s and brings large 
inflows into the centre of the halo. The radial veloc¬ 
ity temporarily becomes positive in a few places due 
to the presence of clumps, and the positive sign then 
describes the flow towards the other clump. This is 
particularly prominent in halo 2 for Z/Zq = 10“"* and 
causes a decline in the mass inflow rate. Repeated 
peaks in the radial velocity indicate the onset of gravi¬ 
tational instabilities which help in the transfer of angu¬ 
lar momentum. The enclosed gas mass in the central 
30 pc of the halo is a few times 10® Mq. The mass 
profile increases with ~ R^ in the central 100 AU and 
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at larger scales linearly increases with radius. Some 
differences in the mass profile are observed between 
metal free and metal poor cases. The enclosed mass is 
a factor of a few lower for Z/Zq = 10 in compari¬ 
son with the metal-free case and is probably a conse¬ 
quence of differences in the thermal structure as well 
as enhanced rotation for Z/Zq = 10 cases. 

Overall, the mass inflow rates at the end of our sim¬ 
ulations seem to be sufficient for forming a supermas- 
sive star at-least for Z/Zq < 10“^ while for higher 
metallicities enhanced rotation as well as fragmenta¬ 
tion may limit accretion onto the central object. How 
long such accretion rates can be maintained depends 
on whether efficient fragmentation takes place or not. 
In the case of efficient fragmentation radiation feed¬ 
back from in situ star formation may shut down ac¬ 
cretion onto the central star and a supermassive star 
cannot form. Even in the absence of fragmentation, an 
enhanced rotation can limit mass accretion as found 
for Z/Zq < 10 in both halos. We discuss these pos¬ 
sibilities in the section below. 


3.2.3. Fragmentation 


To quanti fy the fragment ation, we use the YT 
clump finder (iTurk et al.l201 ih to identify clumps. Our 
clump finding algorithm locates topologically discon¬ 
nected structures. The clumps are considered to be 
bound if the sum of kinetic a nd thermal energy is less 
than the potential energy, see lSmith et alJ ( 20091) . The 
density structure within the central 1 pc of the halo 
is shown in the figure |5] for all cases. There is not 
much substructure found at this scale as the collapse 
is almost isothermal and dust cooling starts to become 
important around a few 1000 AU. Zooming into the 
central 4000 AU (see figure [b]), significant differences 
in the morphology of the halos are obseryed. Particu¬ 
larly, for Z/Zq = 10 the structure is more elongated 
and collapses into a filament due to the more efficient 
cooling compared with other cases where central gas 
clouds are more spherical. ForZ/Z© = 10 two well 
separated grayitationally bound clumps of about a so¬ 
lar mass are formed in halo 2 while in halo 1 the central 
clump is of similar mass, grayitationally bound but the 
second clump is grayitationally unbound and has a 
sub-solar mass. In contrast to this, the clouds are more 
spherical (although not completely) for Z/Zq = 10“^, 
yery similar to a zero metal case. 

For Z/Zq = 10“^, there is more substructure in 
the central part of halo 1 with the possibility of multi¬ 


ple clump formation while halo 2 looks similar to the 
metal free case. Oyerall the clouds are more elongated 
compared to the Z/Z© = 10“^. We also ran a clump 
finder for these cases and found that all clumps in these 
cases are grayitationally unbound and mostly with sub¬ 
solar masses as shown in figure 0 The figure [T] shows 
that only small pockets of gas are cooled by dust cool¬ 
ing. The Z/Zq = 10“® case looks yery similar to the 
metal free case, the temperature of the dense clumps 
of gas is about a thousand K for Z/Zq = 10“^ and it 
declines down to a few hundred K for Z/Zq - 10“^. 
Similar differences in the density structu re of a gas 
cloud are obseryed by IPeters et al.l (2012) employing 
an equation of state with poly tropic index (y <) 1. 

We estimated the maximum scale of thermal insta- 
bi lify (i.e. l,h - CsFnni) fo llowing the criterion giyen 


Inoue & Omukail ( 2015 ) and found that it is com¬ 


parable to the Jeans length. It shows that the res¬ 
olution in our simulations is enough to resolye the 
thermal instability. Dust cooling triggers the Jeans 
instability which results in the formation of multi¬ 
ple clumps. In a self-grayitating cloud, turbulence 
can also induce fragmentation by locally compress- 
ing the gas as obse ryed in an isothermal collapse, see 
Fatif et al. ( 2013ah . Oyerall, central clumps forming 
for Z/Zq = 10 are grayitationally bound and for 
low metallicities all the clumps are grayitationally un¬ 
bound. We caution the reader that grayitationally un¬ 
bound clumps at the initial stages of the collapse do 
not reflect the true amount of fragmentation. To fur¬ 
ther estimate the fate of these clumps, we followed 
the collapse of halo 2 with Z/Zq = 10“^ to densi¬ 
ties of 10“^ g/cm^ by employing 5 additional refine¬ 
ment leyels. In this case, clumps eyentually collapse 
into the centre of the halo. Howeyer, the possibility of 
fragmentation at later stages of the collapse cannot be 
ruled out and it may yary from halo to halo. Eyoly- 
ing these simulations further in time becomes compu¬ 
tationally extremely expensiye due to the large dynam¬ 
ical range coyered in these simulations which results 
in yery short time steps. In the following section, we 
discuss the possible implications of efficient fragmen¬ 
tation if it occurs at later times. 


3.3. Fraction of metal enriched halos 

We estimate the fraction of metal polluted halos 
with Z/Zq < 1 0 from cosmologica l hydrodynam- 
ical simulations Habouzit et al. ( 20161) . The compu¬ 
tational yolume has a size of 10 Mpc and simula¬ 
tions include recipes for star formation and super- 
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nova feedback , a det ailed discussion is provided in 
Habouzit et al. ( 2016h . We identified the halos with 
mass range b etween 2 x 10^ - 10 ^ Mq using the Halo- 
Maker code (ITweed et alJ l2009h at least resolved by 
100 hundred DM particles and computed their mean 
metallicity. Our estimates for the fraction of halos 
with Z/Zq < 10 are shown in figure |9] The num¬ 
ber of metal polluted halos increases with decreasing 
redshift as expected and the fraction of halos polluted 
with Z/Zq < 10“^ is about a factor of 1.5 higher in 
comparison with metal free halos. 


The halos polluted by a trace amount of metals can 
be irradiated by a strong LW flux if they form in the 
close vicinity of a star forming galaxy. Our simula¬ 
tions do not self-consistently take into account the lo¬ 
cal variations in the LW flux required for the direct 
collapse sites and the simulations volume is not large 
enough to capture rare sources which provide a strong 
LW flux of 10^ in terms of J 21 . Assuming that the 
metallicity distribution is however uncorrelated with 
fluctuations in the LW background, if a direct collapse 
is feasible for metallicities up to 10 “^, the resulting 
abundance of the halos is increased by a factor of 1.5 
compared to the purely primordial case. 


4. Implications for black holes formation 


Stellar evolution calculations suggest that the main 
requirement for the formation of direct collapse black 
holes are mass inflow rates of > 0.1 Mo/yr. For 
such accretion rates, stellar evolution differs from nor¬ 
mal stars, as the stellar radius increases with mass 
and no strong UV feedback is produced by the pro¬ 
tostar. Under these conditions, a supermassive star 
of about 10 ^ M 0 can be formed which may later col- 
lapse into a black hole while retaining most of its mas s 


( Hosokawa et al.ll2012l 2013 : Schleicher et al.lEoi3 ). 


In this study, we explored the impact of dust cool¬ 
ing occurring at high densities in the presence of a 
strong LW flux for two halos of a few times 10^ M 0 
and metallicities of Z/Z 0 = 10“^ - 10“^. We resolved 
the collapse down to sub-AU scales and densities of 
about a few times 10'" g/cm^. The presence of large 
mass inflow rates of > 0.1 M 0 /yr suggests that trace 
amounts of metals and dust do not inhibit the forma¬ 
tion of massive objects. At least for Z/Z 0 < 10 '^ the 
mass accretion rates are comparable with an atomic 
cooling case and no strong fragmentation is observed. 
It should be noted that the conditions found here are 
different from minihalos cooled by dust cooling in the 


absence of a strong LW background. Both mass in¬ 
flow rates and available gas mass are about two orders 
of magnitude higher compared to minihalos. Due to 
the extremely large spatial range covered in our simu¬ 
lations, the time step becomes very short and does not 
allow us to evolve simulations long enough to assess 
fragmentation and inflow rates at later times. In the 
following, we consider two possible outcomes assum¬ 
ing that efficient fragmentation occurs at later times. 


In the first case, we consider that efficient frag¬ 
mentation occurs at later times but clumps move in¬ 
ward due to the short mig ration time scale and merge 
with the central clump d lnavoshi & HaimanI 20JA ; 
Latif & Schleicher l20L5albl: ISchleicher et al.l l20L5h . 


In the presence of the large inflow rates of 0.1 Mo/yr 
found here, viscous heating may become important 
in the interior of the disk, it will heat the gas up to a 
few thousand K. Therefore, viscous heating may sup¬ 
press fragmentation by increasing the thermal Jeans 
mass and may help in the fo r mation of massive object s 
( Latif & SchleichCT 2015b : Schleicher et al. 2015 ). 


Moreover, for such large inflow rates the impact of 
UV feedback from stars is expected to be quite weak. 
Even if some of the clumps survive or get ejected via 
3-body processes and form low mass stars, a mas¬ 
sive central object is still expected to form due to the 
clumpy mode of accretion and the stabilisation of the 
disk due to the viscous heating. We argue that in such 
a scenario, massive seed black holes of ~ 10 "^ M 0 can 
form for Z/Z 0 < 10 '^. 


For the second case, we assume that efficient frag¬ 
mentation takes place due to the dust cooling at high 
densities. Then the clumps are unable to migrate in¬ 
ward and a stellar cluster forms. A dense stellar cluster 
is expected to form in such a scenario as dust cool¬ 
ing becomes effective only in the central few 100 - 
1000 AU surrounded by hot gas with a temperature 
around 8000 K. The minimum Jeans mass at high den¬ 
sities where dust cooling operates in our simulations 
is about > 0.1 M 0 and the enclosed mass within a 
radius of a few hundred AU is about 1000 M 0 . As 
a consequence of gravitational collapse, a higher gas 
mass will collapse in the centre over time and may 
reach up to 10 ® M 0 which is about few % of the 
halo mass as often observed in numer ical simulations 
( .Johnson et al. 2009; Latif et al. 2011). For a star for¬ 
mation efficiency of about 10 %, a stellar cluster of 
about 10® M 0 may form. The expected radius of such 
a stellar cluster is about 1000 AU. Various studies 
suggest that such dense clusters may collapse into a 
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massive black hole either via relativistic instabilities 


or stellar dynamical processes (Baumsarte & Shaoiro 

1999; Porteeies Zwart et al 

1999HOmukai et al.ll2008l; 

Deyecchi & Volonteril2009 

Deyecchi et al.l2012h. We 


estimate the mass of a seed black hole forming from 
the core collapse of a cluster in t he following way 
( Portegies Zwart & McMillanll2002 ). 


Mbh = m* + 4 X 10 ^fcMcojliiAc , 


(3) 


where m, is the mass of a massive star in the cluster 
which initiates runaway growth, is the effective frac¬ 
tion of binaries formed dynamically, M^o is the mass 
of the cluster at its birth, InAc is the Coulomb loga- 
rithm and y ~ 1, is the ratio of vario us time scales, see 
Portegies Zwart & McMillanI ( 2002 ). For a star cluster 
of Mco = 10^ Mo, m, = 100 Mq, fc = 0.2, InAc = 10, 
we get Mbh ~ 180 Mq. For a Mco = 10^ Mq, expected 
black holes mass is ~ 900 Mq. Therefore, depending 
on the initial mass of a cluster and the central star, mas¬ 
sive black holes seeds of up to a thousand solar masses 
can be formed in metal-poor halos illuminated by the 
strong LW flux. 


5. Discussion and conclusions 


We have performed high resolution cosmological 
simulations to study the impact of trace amounts of 
metals and dust cooling in massive primordial halos 
irradiated by the strong LW flux. To accomplish this 
goal, we selected two massive halos of a few times 
10 ’ Mq at z > 10 with metallicities of Z/Zq = 10 ““^ - 
10~® and turned on a strong LW radiation of strength 
10^ in terms of J 21 . Our simulations cover a large dy¬ 
namical range by resolving the collapse starting from 
cosmological scales down to scales of sub-AU. To take 
into account the effect of dust and metal line cool¬ 
ing, we extended our previous chemical model and in¬ 
cluded H 2 formation on dust grains, dust grain cooling 
and heating, CII and OI metal lines cooling following 


Omukai et al.l (l2008h . 


Our results show that even in the presence of a 
trace amount of metals and dust, the collapse proceeds 
isothermally with temperatures around 8000 K up to 
densities of about 10“'® g/cm^. Dust cooling becomes 
effective at densities between 10“*® - 10“*’ g/cm’ and 
brings the gas temperature down to 100-1000 K for 
Z/Zq = 10“"* - 10“®, respectively. As expected, dust 
cooling is more efficient and occurs at earlier stages of 
the collapse for higher metallicities in comparison with 
a lower metallicity case. In contrast to the isothermal 


case, H 2 formation on dust grains takes place in the 
core of the halo and its fraction gets boosted. For the 
metal poor cases studied here, metal line cooling is not 
important but could be effective at metallicities higher 
than explored here. 

We found that large inflow rates of 0.1 Mo/yr are 
available for Z/Zq < 10“® while for the higher metal¬ 
licity case the inflow rates start to decline earlier. The 
decline in accretion rate for higher metallicity cases 
is due to the lower gas temperature in the presence 
of dust cooling as well as because of the fragmenta¬ 
tion. No gravitationally bound clumps are found for 
Z/Zq < 10“® by the end of our simulations while for 
Z/Zq = 10“"* clumps are gravitationally bound and 
more massive. The morphology of the central gas 
cloud is significantly different in Z/Zq = 10“"*, as the 
cloud is more elongated compared to other cas es prob¬ 
ably d ue t o the mo r e effic ient dust cooling. LarsonI 
( 1969 ) and PenstonI ( 1969 ) found a self- similar solu- 
tion under isothermal conditions and later Yahil ( 1983h 
showed that it changes with the equation of state by in¬ 
cluding a dependence of the collapse time on the equa¬ 
tion of state parameter y - 1 + The latter sug¬ 

gests that the collapse dynamics will be changed by 
adding dust cooling. Overall, no strong fragmenta¬ 
tion is observed for Z/Zq < 10“® cases although the 
possibility of fragmentation at later stages of the col¬ 
lapse cannot be ruled out. Due to the computational 
constraints, we were unable to evolve simulations for 
longer times and therefore cannot make strong state¬ 
ments about fragmentation occurring at later times. 

The presence of large inflow rates suggests that a 
massive central star is expected to form for Z/Zq < 
10 “® where dust cooling is confined to the cen¬ 
tral 10-100 AU depending on the metallicity. Even 
if the cloud fragments at later stages of the col¬ 
lapse, the clumps are expected to migrate inward and 
merge in the centre as proposed by theoretical work s 


( Inavoshi & Haimanl2014i: Latif & Schleicheil2015b ). 
Viscous heating which becomes particularly important 
in the presence of large inflow rates and rapid rota¬ 
tion may further help in suppressing the fragmentation 
by heating the gas and ev a porating the dus t grain s 
( Latif & Schleicheil 2015b : Schleicher et al.l 2015 ). 
For Z/Zq = 10“"*, efficient fragmentation is expected 
to occur which may lead to the formation of a stellar 
cluster. 

The strength of the LW flux required to keep the 
collapse completely isot hermal in metal-f ree halos is 
about a few times 10"* ( Latif et ^ 2015h and in the 
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present work we employed even a higher value of 
10^ in units of J 21 . The expected number density of 
direct collapse black holes forming in primordial halos 
for Jcrit > 1 O'* is a few orders of magnitude lower than 
the ob served abundance of quasars at z=6 (i.e. 1 per 
Gdc^ dWillott et al~ 201 dl l extrapolating the results of 
Habouzit et alJ ( 2015 ). Even relaxing the constraint of 


primordial gas and allowing the direct collapse black 
holes to form in halos with Z/Z© < 10“^ in the pres¬ 
ence of a strong LW flux only increases their abun¬ 
dance by a factor of two. Even under these conditions, 
the formation sites for the direct collapse black holes 
remain rare. 

Eor metallicities higher than explored here, cooling 
due to the fine structure lines of CII and OI becomes 
important at low densities about 10^“° g/cm^ even in 
the presence of a strong LW flux and cools the gas 
down to a few 10 K. In such cases, the thermal insta¬ 
bility induces fragmentation and det ermines the mass 
spectr um of the resulting star cluster (llnoue & Omukai 


20L5h . 
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Fig. 1.— Figure shows the temperature, H 2 , Ft' and e' fractions plotted against the density from one-zone test. The 
green, blue and red lines represent Z/Z© = 10'®, 10 ® and 10'^, respectively. 
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Fig. 2.— Density-temperature phase diagram for Z/Zq = 10“^, 10'^ and 10 from top to bottom, respectively. The 
left panel represents halo 1 and the right panel halo 2. The gas is initially heated up to a few times 10"^ K, collapses 
isothermally up to the densities of 10 g/cm^ in the presence of a strong LW flux and the dust cooling becomes 
effective at > 10 g/cm^ densities depending on the amount of metallicity. Colour-bars show the amount of gas in 
solar masses for a given density. 
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Radius [AU] Radius [AU] 

Fig. 3.— Radially averaged and spherically binned profiles of various quantities for halo 1 (HI) and halo 2 (H2) 
are shown here. The green, blue and red lines represent Z/Z© = 10 10“^ and 10 respectively. The solid lines 

represent halo 1 and the dashed lines show halo 2. The vertical black dashed line corresponds to the Jeans length. An 
isothermal case for halo 1 with zero metallicity is also plotted as a reference and is represented by the black solid line. 
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Fig. 4.— The ratio of rotational to Keplerian velocity for halo 1 (HI) and halo 2 (H2) is shown here. The green, blue 
and red lines represent Z/Z© = 10 10“^ and 10“®, respectively. The solid lines represent halo 1 and the dashed lines 

show halo 2. An isothermal case for halo 1 with zero metallicity is also plotted as a reference and is represented by 
the black solid line. The vertical black dashed line corresponds to the Jeans length. 
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Fig. 5.— Average gas density along x-axis for the central 1 pc of a halo. Each row represents a halo (halo 1(H1) on 
top and halo 2 (H2) on bottom) and each column represent metallicity (increasing from left to right). HlZe-6, HlZe-5 
and HlZe-4 represent Z/Z© = 10“^, 10“^ and 10“^ for halo 1, respectively. Similarly, H2Ze-6, H2Ze-5 and H2Ze-4 
represent Z/Zq = 10 10“^ and 10 for halo 2. 
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Fig. 6.— Average gas density along x-axis for the central 4000 AU of a halo. Each row represents a halo (halo 1 on 
top and halo 2 on bottom) and each column represent metallicity (increasing from left to right). 
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Fig. 7.— Average temperature weighted by the gas density along x-axis for the central 4000 AU of a halo. Each row 
represents a halo (halo 1 on top and halo 2 on bottom) and each column represent metallicity (increasing from left to 
right). 
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Fig. 8.— The ratio of clump mass to the Jeans mass is plotted against the clump masses for various metallicities. The 
left panel halo 1 and the right panel halo 2. The red symbols represent the clumps forming for Z/Z© = 10~®, the blue 
for Z/Zq = 10“^ and the green for Z/Zq = 10“^. The clumps for Z/Zq < 10“^ are gravitationally unbound while the 
central clumps for Z/Zq = 10“^ cases are gravitationally bound and more massive. In general, the ratio of clump mass 
to the Jeans mass increases with clump masses. 
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Fig. 9.— Fraction of halos with metallicities below a given value (see legend) as a function of redshift. In this figure, 
we only show the fraction of halos with masses between 2 x 10^ - 10^ Mq in a computational box of size 10 Mpc for 
Z/Zq < 10“^. The fraction of halos with Z/Zq < 10“^ is about 1.5 times higher than the metal-free halos. 
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